function EPSC_train()
% based on the EPSC2 data for train responses, compute the mean and stdev
% of the EPSC amplitudes through the train. 4-6/2008.
%


global CONTROL CURRENT DFILE
sf_list = getmainselection;
use_d = 0;
for sf = sf_list

    [DFILE, err] = analysis_setup(DFILE, sf); % get data and parameters.

    if(err ~= 0)
        return;
    end;
    bl = number_arg(CONTROL(sf).durho); % get the baseline from the holding duration
    st = CONTROL(sf).stim_time;
    if(ischar(st))
        st = number_arg(st); %#ok<NASGU>
    end;
    tm = make_time(DFILE);
    dsr = 0.001*(DFILE.rate(1)*DFILE.nr_channel(1));
    if(isfield(CONTROL, 'EPSC2'))
        d = CONTROL(sf).EPSC2;
    else
        d=[];
    end;
    if(~isempty(d)  && use_d)
        t0 = find(tm(1,:) >= d.tstart, 1);
        t1 = find(tm(1,:) >= d.tend, 1);
    else
        t0 = 1;
        t1 = size(tm, 2);
    end;
    pro = CONTROL(sf).protocol;

    par = CONTROL(sf).Parameters;

    r = make_struct(par, 'delay 5, np 20, freq 20, tau_init 40, nrep 4');
    freq = r.freq;
    tau_init = r.tau_init;
    nrep = r.nrep;
    np = r.np;
    delay = r.delay;

    

    acurr = CURRENT;
    da = diff(CURRENT, 1, 2);
    ua = reshape(da, numel(da), 1);
    umax = max(ua); umin = min(ua);
    uh = umin:100:umax;
    dh = hist(ua, uh); % convert to histogram
    duh = find(uh > max(uh)/2);
    thr = min(uh(duh(1)));
    thr = 4000;
    [acurr, st] = FP_artsupp(acurr, DFILE, sf, 1, thr); % suppress the artifacts if stimulating.....

    
    tisi = 1000/freq;
    st = delay - tisi + cumsum(tisi*ones(1, np));
    rdelay = max(st);
    if(strcmp(DFILE.filename, '08sep08b') && DFILE.frec == 352)
        sequence ='6;2100/15l';
    else
        sequence = '5;2000/14l';
    end;
    strec = seqparse(sequence);
    strec = strec + rdelay;
    [acurr, str] = FP_artsupp(acurr, DFILE, sf, 1, thr, strec);
    
    %     thr = 18000;
    %     while(isempty(st)) % do this until we have some artifacts suppressed - only a guess...
    %         thr = thr/10;
    %         [acurr, st] = FP_artsupp(acurr, DFILE, sf, 1, thr); % suppress the artifacts if stimulating.....
    %     end;
    if(isempty(freq) || isnan(freq)) % we wait unitl here so we have a valid st to work with
        freq = 1000/mean_var(diff(st(1:end-1)));
    end;

    basel = 1:(floor(bl/dsr+0.5)-1); % subtract the baseline FIRST
    for i = 1:size(acurr, 1)
        acurr(i,:) = acurr(i,:) - mean_var(acurr(i, basel));
    end;
    tisi = 1000/freq; % mean(diff(st));
    %   stdisi = std(diff(st));
    [bad, ibad] = find(diff(st) > 1.1*tisi);
    if(~isempty(ibad))
        st(ibad+1:end) = NaN;
        st=st(find(~isnan(st)));  %#ok<FNDSB>
    end;
    %    tisi = 1/freq; % mean(diff(st)); % re-evaluate mean interval..
    npulse = length(st);
    if(~isempty(d) && use_d)
        n = length(d.eventtime);
    else
        n = size(tm, 1);
    end;
    amps = zeros(n, npulse);
    slowamp = zeros(n, npulse);
    tn = amps;
    %t = tisi*(1:npulse);
    thisfigure = newfigure('EPSC_Train', 'EPSC Train Analysis');
    set(gcf, 'color', 'white');
    subplot('Position', [0.07, 0.6, 0.8, 0.37]);

    % acurr = mean(acurr(:,t0:t1));
    plot(tm(1,t0:t1), acurr);
    xlabel('T (ms)');
    ylabel('I (pA)');
    xmax = min(max(tm(1,t0:t1)), max(strec)+100);
    set(gca, 'Xlim', [0 xmax]);
    set(gca, 'box', 'off');
    if(~isempty(d) && use_d)
        for i = 1:n
            j = floor(d.eventtime{i}/tisi);
            if(min(j) == 0)
                j = j + 1;
            end;
            amps(i,j) = d.eventamp{i};
            tn(i,j) = j*tisi;
        end;
    else
        kts = floor(st/dsr); % get start times
        kstep = floor(tisi/dsr);  % time between (for measurement window)
        for i = 1:size(acurr, 1)
            for j = 1:npulse
                slowamp(i,j) = -mean_var(acurr(i, kts(j)-1:kts(j)));
                ampe = kts(j)+kstep-1;
                if(ampe > size(acurr, 2))
                    ampe = size(acurr, 2) - 1;
                end;
                amps(i,j) = -(min(acurr(i, kts(j): ampe)+slowamp(i,j)));
            end;
        end;
    end;
    fsum = slowamp./(amps+slowamp); % fraction of current due to summation
    %
    % do the same for the recovery
    ktsr = floor(strec/dsr); % get start times
    kstep = floor(tisi/dsr);  % time between (for measurement window) - keep same for recovery
    for i = 1:size(acurr, 1)
        for j = 1:length(strec)
            slowamprec(i,j) = -mean_var(acurr(i, ktsr(j)-1:ktsr(j)));
            ampe = ktsr(j)+kstep-1;
            if(ampe > size(acurr, 2))
                ampe = size(acurr, 2) - 1;
            end;
            ampsrec(i,j) = -(min(acurr(i, ktsr(j): ampe)+slowamprec(i,j)));
            %                fprintf(1, 'i=%d, j=%d: (t=%6.1f) I = %10.3f\n', i, j, tm(i, ktsr(j)), ampsrec(i,j));
        end;
    end;





    subplot('Position', [0.07, 0.3, 0.4, 0.24]);
    plot(st, amps, 'k.')
    xlabel('T (ms)');
    ylabel('I_{EPSC} (pA)');
    upsc = get(gca, 'YLim');
    set(gca, 'YLim', [0 upsc(2)]);
    set(gca, 'box', 'off');
    hold on;
    drawnow
    npulse = length(st);
    eavg = zeros(1, npulse);
    estd = zeros(1, npulse);
    seavg = eavg;
    sestd = estd;
    fsumave = eavg;
    fsumstd = eavg;
    for i = 1:npulse
        [eavg(i) estd(i)] = mean_var(amps(:,i));
        estd(i) = sqrt(estd(i)/length(amps(:,i)));
        [seavg(i) sestd(i)] = mean_var(slowamp(:,i));
        sestd(i) = sqrt(sestd(i)/length(slowamp(:,i)));
        [fsumave(i) fsumstd(i)] = mean_var(fsum(:,i));
        fsumstd(i) = sqrt(fsumstd(i)/length(fsumstd(:,i)));
    end;
    %length(st)
    %length(eavg)

    startfitpulse = 2;
    stf = (st-st(1)); % zero starts at first pulse (simplifies interpretation)
    tf = (stf(startfitpulse):0.1:stf(end)); % and the result array starts at the first pulse

    model=50;


    % pars:    DC     A0   Tau1
    pmask =   [1,       1,         1];
    lbound =  [-25000, -25000,     1];
    ubound =  [ 25000,  25000,   500];
    lam(1) = 0;
    lam(2) = 1000;
    lam(3) = tau_init;
    order=length(lam)-1;
    alpha = 0;
    beta = 0;
    maxiter = 2000;
    %
    [c,lam, goodfit]=curve_fitting(stf(startfitpulse:end)-stf(1), eavg(startfitpulse:end), ...
        'levenberg','cubic', model, order, lam, pmask, lbound, ubound, alpha, beta, maxiter); %#ok<NASGU>

    [f, yf]=fit_func(lam, tf, zeros(size(tf)), model, alpha, 0); % calculate fit
    my_errorbar(st, eavg, -estd, estd,'ro', 1);
    plot(tf+st(1), yf, 'g-');
    u = get(gca, 'Xlim');
    if(u(1) < 0) % keep axis scaling sane.
        u(1) = -10;
        set(gca, 'Xlim', u);
    end;

    % subplot for recovery
    subplot('Position', [0.54, 0.30, 0.4, 0.24]);
    for i = 1:size(acurr,1) % for all the records
        j = floor((i-1)/nrep)+1; % repeat counter index
        n = mod(i, nrep) + 1;
        irecover(j,n) = ampsrec(i, j);
 %       fprintf(1, 'i: %d j: %d  irec: %8.1f\n', i, j, irecover(i));
    end;
    recover_delay = strec - rdelay;
    hrec = semilogx(recover_delay, irecover);
    hold on;
    set(hrec, 'Marker', 's', 'MarkerFaceColor', 'k', 'MarkerEdgeColor', 'k', 'MarkerSize', 3);
    set(hrec, 'LineStyle', 'none');
    plot([recover_delay(1),recover_delay(end)], [eavg(1), eavg(1)], 'r-')

    set(gca, 'Ylim', [0 upsc(2)]);
    set(gca, 'Xlim', [min(recover_delay)/2, max(recover_delay)*1.1]);
    
    subplot('Position', [0.07, 0.07, 0.4, 0.15]);
    my_errorbar(st, seavg, -sestd, sestd,'ks', 1);
    xlabel('T (ms)');
    ylabel('I_{slow} (pA)');
    set(gca, 'box', 'off');
    % use right axis for separate plot on this same graph
    ax1 = gca; % get the first axis
    ax2 = axes('Position',get(ax1,'Position'),...
        'XAxisLocation','top',...
        'YAxisLocation','right',...
        'Color','none',...
        'XColor','b','YColor','b', ...
        'XTickLabel', '');
    hold on;
    he = my_errorbar(st, fsumave, -fsumstd, fsumstd, 'ob', 1);
    set(ax2, 'Ylim', [-0.1 1.25]);
    ylabel('Fraction (slow)');
    set(he, 'Parent', ax2);
    u = get(gca, 'Xlim');
    if(u(1) < 0) % keep axis scaling sane.
        u(1) = -10;
        set([ax1 ax2], 'Xlim', u);
    end;
    subplot('Position',[0.65,0.0,0.32,0.25])
    axis([0,1,0,1])
    axis('off')
    d.channel = 1;
    text(0,0.9,sprintf('%-12s R[%d:%d] Channel: %d', ...
        DFILE.filename, DFILE.frec, DFILE.lrec,  d.channel), 'Fontsize', 11);
    text(0, 0.8, sprintf('%-25s', CONTROL(sf).protocol), 'Fontsize', 11);
    text(0, 0.7, sprintf('Fit: A0 = %7.3f   A1 = %7.3f ', lam(1), lam(2)), 'Fontsize', 11);
    text(0, 0.6, sprintf(' Tau = %7.3f', lam(3)), 'Fontsize', 11);
    text(0, 0.5, sprintf(' Cell Type: %s', CONTROL(sf).CellType), 'Fontsize', 11, 'FontWeight', 'bold');
    text(1, 0.1, sprintf('EPSC_Train V1.2 6/03/08 pbm'), 'Fontsize', 9, 'Interpreter', 'none',...
        'HorizontalAlignment', 'right');

    if(~isempty(freq))
        EPSC_Train.Freq = freq; % rounds off to nearest 10 Hz... cheat
    else
        EPSC_Train.Freq = 0;
    end;
    EPSC_Train.amps = amps;
    EPSC_Train.slowamp = slowamp;
    EPSC_Train.stim = st;
    EPSC_Train.fit_y  = yf;
    EPSC_Train.fit_x = tf;
    EPSC_Train.fitpar = lam;

    EPSC_Trecover.Freq = freq;
    EPSC_Trecover.seq = sequence;
    EPSC_Trecover.Tdelays = strec;
    EPSC_Trecover.rdelay = rdelay;
    EPSC_Trecover.Iamps = irecover;

    CONTROL(sf).EPSC_Train = EPSC_Train;
    CONTROL(sf).EPSC_Recovery = EPSC_Trecover;
    [p f] = fileparts(DFILE.filename);
    figurefile = sprintf('%s/%s_R%d_%dHz.pdf', DFILE.path, f, DFILE.frec, EPSC_Train.Freq);
    saveas(thisfigure, figurefile, 'pdf');
    fprintf(1, 'EPSC_Train: plot saved to %s\n', figurefile);
end;

